namespace Eigen {

namespace internal {

    template <typename Scalar>
    void r1updt(Matrix<Scalar, Dynamic, Dynamic>& s,
                const Matrix<Scalar, Dynamic, 1>& u,
                std::vector<JacobiRotation<Scalar>>& v_givens,
                std::vector<JacobiRotation<Scalar>>& w_givens,
                Matrix<Scalar, Dynamic, 1>& v,
                Matrix<Scalar, Dynamic, 1>& w,
                bool* sing)
    {
        typedef DenseIndex Index;
        const JacobiRotation<Scalar> IdentityRotation = JacobiRotation<Scalar>(1, 0);

        /* Local variables */
        const Index m = s.rows();
        const Index n = s.cols();
        Index i, j = 1;
        Scalar temp;
        JacobiRotation<Scalar> givens;

        // r1updt had a broader usecase, but we don't use it here. And, more
        // importantly, we can not test it.
        eigen_assert(m == n);
        eigen_assert(u.size() == m);
        eigen_assert(v.size() == n);
        eigen_assert(w.size() == n);

        /* move the nontrivial part of the last column of s into w. */
        w[n - 1] = s(n - 1, n - 1);

        /* rotate the vector v into a multiple of the n-th unit vector */
        /* in such a way that a spike is introduced into w. */
        for (j = n - 2; j >= 0; --j)
        {
            w[j] = 0.;
            if (v[j] != 0.)
            {
                /* determine a givens rotation which eliminates the */
                /* j-th element of v. */
                givens.makeGivens(-v[n - 1], v[j]);

                /* apply the transformation to v and store the information */
                /* necessary to recover the givens rotation. */
                v[n - 1] = givens.s() * v[j] + givens.c() * v[n - 1];
                v_givens[j] = givens;

                /* apply the transformation to s and extend the spike in w. */
                for (i = j; i < m; ++i)
                {
                    temp = givens.c() * s(j, i) - givens.s() * w[i];
                    w[i] = givens.s() * s(j, i) + givens.c() * w[i];
                    s(j, i) = temp;
                }
            }
            else
                v_givens[j] = IdentityRotation;
        }

        /* add the spike from the rank 1 update to w. */
        w += v[n - 1] * u;

        /* eliminate the spike. */
        *sing = false;
        for (j = 0; j < n - 1; ++j)
        {
            if (w[j] != 0.)
            {
                /* determine a givens rotation which eliminates the */
                /* j-th element of the spike. */
                givens.makeGivens(-s(j, j), w[j]);

                /* apply the transformation to s and reduce the spike in w. */
                for (i = j; i < m; ++i)
                {
                    temp = givens.c() * s(j, i) + givens.s() * w[i];
                    w[i] = -givens.s() * s(j, i) + givens.c() * w[i];
                    s(j, i) = temp;
                }

                /* store the information necessary to recover the */
                /* givens rotation. */
                w_givens[j] = givens;
            }
            else
                v_givens[j] = IdentityRotation;

            /* test for zero diagonal elements in the output s. */
            if (s(j, j) == 0.)
            {
                *sing = true;
            }
        }
        /* move w back into the last column of the output s. */
        s(n - 1, n - 1) = w[n - 1];

        if (s(j, j) == 0.)
        {
            *sing = true;
        }
        return;
    }

}  // end namespace internal

}  // end namespace Eigen
